# generate graphs of coefficients #

  # estimate effects (raw) #

    # run regressions #

      coef = rep(0,length(depvars)) # vector to store coefficients #

      se = rep(0,length(depvars)) # vector to store standard errors #

    # run regressions #

      for (i in 1:length(depvars)) {
  
        f = as.formula(paste(depvars[i],treatment, sep = " ~ "))
        model = lm_robust(f, data = agro, se_type = "HC1")
        coef[i] = model$coefficients[2]
        se[i] = sqrt(diag(vcov(model))[2])
  
      }

    # store coefficients and standard errors as text #

      text = rep("",length(depvars)) # vector to store text #
      beta='\u03b2'

      for (i in 1:length(text)) {
  
        text[i] = paste0(beta," (se) = ", 
                         formatC(coef[i], format = 'f', flag='0', digits = 2)," (",
                         formatC(se[i], format = 'f', flag='0', digits = 2),")")
  
      }

    # estimate effects (normalized) #

      agro = agro %>% mutate_at(depvars,scale)
      
      # data to store results #
    
        spec =  as.data.frame(rep(3:1, each = length(depvars)))
        number = as.data.frame(c(1:length(depvars)))
        colnames(spec) = "model"
        colnames(number) = "depvar"

      # specification 1 #
      
        coef = rep(0,length(depvars)) # vector to store coefficients #
        se = rep(0,length(depvars)) # vector to store standard errors #
        for (i in 1:length(depvars)) {
          f = as.formula(paste(depvars[i],treatment, sep = " ~ "))
          model = lm_robust(f, data = agro, se_type = "HC1")
          coef[i] = model$coefficients[2]
          se[i] = sqrt(diag(vcov(model))[2])
        }
        results1 = cbind(number,coef,se)
        results1 = results1 %>% mutate(icmax = coef+1.96*se, icmin = coef-1.96*se)
        
      # specification 2 #
        
        coef = rep(0,length(depvars)) # vector to store coefficients #
        se = rep(0,length(depvars)) # vector to store standard errors #
        for (i in 1:length(depvars)) {
          f = as.formula(paste0(depvars[i]," ~ ",treatment,c(" + yob + sex + res + main + coop + area")))
          model = lm_robust(f, data = agro, se_type = "HC1")
          coef[i] = model$coefficients[2]
          se[i] = sqrt(diag(vcov(model))[2])
        }
        results2 = cbind(number,coef,se)
        results2 = results2 %>% mutate(icmax = coef+1.96*se, icmin = coef-1.96*se)
        
      # specification 3 #
        
        coef = rep(0,length(depvars)) # vector to store coefficients #
        se = rep(0,length(depvars)) # vector to store standard errors #
        for (i in 1:length(depvars)) {
          f = as.formula(paste0(depvars[i]," ~ ", treatment, c(" + yob + sex + res + main + coop + area")," + ",lagvars[i]))
          model = lm_robust(f, data = agro, se_type = "HC1")
          coef[i] = model$coefficients[2]
          se[i] = sqrt(diag(vcov(model))[2])
        }
        results3 = cbind(number,coef,se)
        results3 = results3 %>% mutate(icmax = coef+1.96*se, icmin = coef-1.96*se)
        
      # data #
      
        results = rbind(results1,results2,results3)
        results = cbind(spec,results)
        results$text = text

  # plot coefficients #
        
    graph = ggplot(results, aes(x = depvar, y = coef, colour = as.factor(model))) + 
      geom_hline(yintercept = 0, colour = "#323232", lty = 2) + 
      geom_point(size = 2, position = position_dodge(width = 1/2)) +
      geom_linerange(aes(ymin = icmin, ymax = icmax), size = 1, position = position_dodge(width = 1/2)) +
      theme_linedraw() + 
      xlab("") + ylab("Standardized Effects") + 
      theme(plot.title = element_text(size=14, family="sans", hjust = 0.5, face = "bold"),
            axis.title.x = element_text(size=13, family="sans", face = "bold", colour = "#323232"), 
            axis.title.y = element_text(size=13, family="sans", face = "bold", colour = "#323232"),
            axis.text.x = element_text(size=11, family="sans", colour = "#323232"), 
            axis.text.y = element_text(size=12, family="sans", colour = "#323232"),
            #panel.grid.major =  element_blank(), 
            panel.grid.minor = element_blank(),
            axis.line = element_line(colour = "#323232", size = 0.6), 
            axis.ticks = element_line(colour = "#323232", size = 1),
            plot.margin = margin(1, 0.5, 0.5, 0.5, "cm"),
            legend.position = "bottom", legend.title = element_blank(), 
            legend.text = element_text(size=12)
          ) +
          scale_x_continuous(breaks=1:8, labels = c("log(Expenditures)",
                                                    "Pesticide Use (0/1)", "Tractor Use (0/1)", 
                                                    "Conserv. Practices (#)", "Manag. Practices (#)", 
                                                    "Rotational Grazing (0/1)", "Pasture Restoration (%)", 
                                                    "Technical Assistance (0/1)")) + 
          scale_y_continuous(breaks = c(-0.4,0,0.4,0.8,1.2,1.6,2.0), limits = c(-0.4,2.4)) + coord_flip() +
          geom_text(aes(x=depvar, y=2, 
                        label = text),
                    size = 4, hjust=0.35, vjust=-0.14, colour = "#323232", family="sans") + 
          scale_color_manual(values=c("#d7191c","#fdae61","#2c7bb6"), 
                             breaks = c("3","2","1"),
                             labels = c("Bivariate","Demographic Controls","Baseline Dep. Var."),
                             name="")
    
    